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Abstract 

Since the discovery of the figure-8 orbit for the three-body problem 
[Moore 1993] a large number of periodic orbits of the n-body prob- 
lem with equal masses and beautiful symmetries have been discovered. 
However, most of those that have appeared in the literature are either 
planar or are obtained from perturbations of planar orbits. Here we 
exhibit a number of new three-dimensional periodic n-body orbits with 
equal masses and cubic symmetry, including some whose moment of 
inertia tensor is a scalar. We found these orbits numerically by mini- 
mizing the action as a function of the trajectories' Fourier coefficients. 
We also give numerical evidence that a planar 3-body orbit first found 
in [Henon, 1976], rediscovered by [Moore 1993], and found to exist for 
different masses by [Nauenberg 2001], is dynamically stable. 

It is a pleasure to dedicate this paper to Philip Holmes. 
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1 Introduction 



Inspired by the work of Holmes and his collaborators on the classification of 
periodic orbits in three-dimensional flows according to their knot and braid 
types [ini [HI [121 [8], in 1993 Moore [15] considered periodic orbits of the n- 
body problem in the plane, and asked which braids can appear in their space- 
time trajectories. We can search for such orbits numerically by starting with 
a fictional orbit with a desired symmetry and topology, and then minimizing 
the action S = J Cdt using gradient descent, where C = — \^ is the 
Lagrangian and K and V are the kinetic and potential energies respectively. 
If r parametrizes the progress of this descent, then with each step 6t we 
increment the position Xj(t) of the ith mass at each time t by 5xj(t), where 

S S 

5x,(t) = -Sr^—^ = 5r(m,x,(t) - F,(t)) (1) 

where Fj(t) = —dV/d:x.i(t) is the force acting on the ith mass at time t, and 
Xi(i) = (Pxi/dt"^ is its acceleration. Each time we update Xj(t) in this way, 
the action 5" decreases. This process can lead to one of three outcomes as r 
increases towards infinity: escape, in which Xj(t) tends to infinity for some 
i and some t; collision, in which Xj(t) = Xj(t) for some i ^ j at some t; or 
success, in which the Xj(t) converge to an orbit which is a (local) minimum 
of the action, and therefore a genuine trajectory of the equations of motion 
where F, = mjXj for all i. 

For homogeneous potentials V{r) oc r", it is easy to show [15J that escape 
cannot occur if a < 2 unless the braid can be separated into two pieces, and 
collision cannot occur (starting with a fictional orbit with finite action) if 
a < —2 since the action of a colliding trajectory in that case is infinite. 
Thus all braid types occur for a < —2, including the so-called "strong force" 
a = —2, but only some subset of the possible braids will exist for, say, 
Newtonian gravity, where a = —1. 

Performing this minimization on discretized orbits, Moore [15] found a 
number of new periodic 3-body orbits in the plane, with a variety of braid 
types and symmetries: in particular, an orbit in which three equal masses 
chase each other around a figure-8 as shown in Figure [T], describing a clas- 
sic three-strand braid in space-time. Surprisingly, simulations indicated that 
this orbit is stable to small perturbations. Six years later, this orbit was re- 
discovered by Chenciner and Montgomery P| , who provided a rigorous proof 
of its existence. More precise numerical work was done by Simo [20] and by 
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Figure 1: The figure-8 orbit for three equal masses. The points show the 
masses at five equal time intervals. 



Nauenberg |17j; Simo identified many new orbits, called "choreographies," 
where masses follow each other along a fixed trajectory, and provided addi- 
tional numerical evidence that the figure-8 orbit is stable to perturbations in 
the plane [2T] . 

However, naively performing gradient descent in the variables Xi{t) can 
introduce numerical problems. In particular, discretizing the second time 
derivative as x(t) = x(t + l) — 2x(t) + x(t — 1) produces a "zig-zag" instability 
in which the orbit jumps back and forth with alternate steps in t and r. 
Moreover, this method can only find local minima of the action, and thus is 
useless if the solution is another type of extremum. 

An alternative approach was considered by Nauenberg [T7] which capital- 
izes on the periodicity of the coordinates Xi{t), by expanding these variables 
in a Fourier series, and then finding extrema of the action as a function 
of the Fourier coefficients of the trajectory. These extrema can be found by 
a gradient procedure which is analogous to the approach in ([T]), but avoids 
the instabilities associated with that method, and is not restricted to finding 
minima of S. For forces associated with the homogeneous potential V oc r°', 
the orbits scale as x{t) = T^/*^^~"^/(t/T) where T is the period, and for our 
discussion we choose T = 2tt. We substitute in S the Fourier expansion 



for each mass, and similarly for y{t) and z{t) (we suppress the subscript i 




(2) 



k 
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and ignore cosine terms for simplicity). Then we have, to first order in a 
small change 6ak in the Fourier coefficients a^, 

= E (3) 

As in ([1]) we perform gradient descent, but now with the rate of descent in 
each coefficient at controlled by a parameter 5Tk which can depend on k: 

Sak = . (4) 

oak 

In particular, substituting K = {1 / 2)m{dx / dtY for the kinetic energy in S 
yields 

^o-k = —^T~k I rnk'^ak + / dt — — sin kt ] . (5) 
V Jo dx J 

Starting with given values of the Fourier coefficients and iterating with 
finite steps 6Tk leads to new values of which decrease the action S, lead- 
ing as before to a local minimum of S. Moreover, by setting some Fourier 
coefficients to zero or imposing equalities between the coefficients of differ- 
ent masses, we can restrict our search to orbits with certain symmetries, 
and by choosing the starting trajectory we can search for orbits with certain 
topologies. 

We recover the naive gradient descent equation of [Uj by setting Sr^ = 6t 
independent of k. Note from ([2]) that 5x(t) = J2k sin kt; then substitut- 
ing ([5]) yields ([1]). However, allowing Srk to depend on k allows us to avoid 
the above-mentioned zig-zag instability. This instability, which is associated 
with discretizing the second time derivative of x, is reflected here in the 
quadratic term in k on the right hand side of (jS]); if we set 6Tk = S/mk"^, 
where 5 is a small positive constant, this removes this high-frequency insta- 
bility, and gives a procedure which converges to a smooth trajectory whose 
Fourier coefficients tend to zero as /c — > oo. In addition to removing the 
high-frequency instability, allowing St^ to depend on k can, in some cases, 
allow us to find solutions which are extrema but not minima of S: if 6Tk is 
negative for some values of k, then we can sometimes find orbits which are 
saddle points of the action. 

Using this approach, Nauenberg [171 EH] found extensions in three dimen- 
sions of the figure-8 which have nonzero angular momentum about the two 
axes of symmetry in the original plane of the eight, and are periodic in a 
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rotating frame; in particular, for the rotating-8 around the y axis, these were 
found with negative values of 5Tj. for k = 1 and A; = 3. Likewise, he found 
extended orbits with angular momentum normal to this plane, which remain 
planar. For the case where the axis of rotation is along the horizontal sym- 
metry axis of the figure-8, a proof for the existence of these orbits was given 
by C. Marchal who showed that these orbits interpolate between the 
figure-8 (which has zero total angular momentum) and the Lagrange orbit 
(where three masses co- rotate around a common center); see also [4]. 

In this paper we present several new three-dimensional periodic orbits. 
These orbits have symmetry that is fully three-dimensional, unlike, for in- 
stance, the "hip-hop" orbit of [2], which is obtained by perturbing a planar 
Lagrange orbit and has only dihedral symmetry. In particular, we find a 
family of orbits with 4m masses where m is an odd integer, with m masses 
on each of four loops which are related by cubic symmetry. These orbits have 
zero total angular momentum. Moreover, if m is a multiple of 3 they have 
cubic symmetry at all times, and their mass distribution is spherically sym- 
metric to second order. Although we have not attempted to do this, we note 
that it should be possible to prove the existence of these orbits rigorously 
using techniques like those in [TU [H [5] . 

We then discuss another orbit we call the "criss-cross," in which two 
masses orbit each other while the third rotates around them both in the op- 
posite direction. In fact, this orbit was first found by Henon as a member 
of a finite angular momemtum family of orbits periodic in rotating frames, 
extended from the zero angular momentum Schubart orbit It was redis- 
covered by Moore [15] by searching for orbits with a particular braid type, 
and later extended to different masses by Nauenberg [T7j. Recently, its ex- 
istence was proved rigorously by Chen [I] as part of a family of retrograde 
orbits. Here we give numerical evidence that this orbit is dynamically sta- 
ble in three dimensions, allowing the tantalizing possibility that it can be 
realized in astronomical systems. 

2 Orbits with cubic symmetry 

One source of inspiration for symmetric orbits in three dimensions consists of 
polyhedra. For example, truncating the corners of a cube until a pair of tri- 
angular faces meet at the midpoint of each edge yields an Archimedean solid, 
the cuboctahedron, shown on the left in Figure [21 Its 24 edges can be thought 
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Figure 2: Left, a cub octahedron, and its four hexagonal "great circles"; right, 
the 12-mass orbit with cubic symmetry. 

of as composing 4 hexagonal "great circles", and these can be oriented in 
such a way that half of the triangular faces are oriented counterclockwise 
and the others are oriented clockwise. It takes only a little imagination to 
suppose that masses can start at the the corners of the cuboctahedron, and 
then travel along these great circle loops in these directions: first orbiting 
around the counterclockwise triangles, and then around the clockwise ones. 

The simplest such orbit has 4 masses, one on each of the four loops. At 
all times, the masses are related to each other by vr-rotations about the x, y 
and z axes: that is, by the transformations 

{x,y,z) {x,y,z),{x,-y,~z),{-x,y,-z),{-x,-y,z) . (6) 

Each mass rotates counterclockwise around a diagonal axis, (+1,+1,+1), 
(+1,-1,-1), (— 1,+1,— 1), or (— 1,— 1,+1). The total angular momentum 
is zero, a fact which we can immediately see from symmetry, since no nonzero 
vector J is fixed under the four transformations of ([6]). 

To find this orbit numerically, we construct circular orbits which satisfy 
these symmetry conditions, and use these as the starting point for the pro- 
cess of action minimization described in the Introduction. For example, the 
coordinates of a circular orbit with unit radius normal to the diagonal (1,1,1) 
are given by 

x{t) = sint, y{t) = sin(t + 27r/3), and z{t) = sin(t + Att/3) . (7) 
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Since our minimization process preserves the symmetries of the orbit, for all 
times t the actual trajectory of this mass has the form 



x{t) = fit), yit) = fit + 27r/3), and ^(t) = fit + 47r/3) . (8) 



and the trajectories of the other masses are related to this one by the trans- 
formations ([6]). 

We now expand fit) in a Fourier series 



and our goal is to find the Fourier coefficients and hk- However, if we set 
/(O) = so that the initial position of this mass is in the y — z plane, by time 
reversal symmetry we have fit) = — /(— t) and the bk are identically zero. 
Moreover, these orbits also have the symmetry that after half a period each 
mass is diametrically opposite to its previous position, so fit + n) = —fit) 
and the Fourier coefficients are zero except for odd k. 

Minimizing the action, starting with these circular orbits and using this 
restricted set of Fourier coefficients, gives the periodic orbit whose Fourier 
coefficients are given in the first column of Table [H To verify adequate 
convergence of our Fourier sums, we compared the accelation of the resulting 
trajectory with the force for a series of values of t. By considering frequencies 
up to A; < 27 we found agreement between the acceleration and the force up 
to 5 decimal figures. Higher accuracy can be achieved by including higher 
frequencies in the Fourier sums. We also integrated the equations of motion 
starting with the initial values we obtained, namely 



X = (0, -0.69548, 0.69548), i = (0.87546, -0.31950, -0.31950) (9) 



using the Runge-Kutta algorithm. We found good agreement during a single 
period; over longer times the integrated orbit rapidly diverges, indicating 
that this periodic orbit is unstable. 

Next, we consider orbits with m masses on each of the four loops. Now, in 
addition to the symmetries described above, these masses are related by the 
Lagrange symmetry, in which the j masses are related by time translation, 
t ^ t + 2jTT /m for j = 0,...,m — 1. 

In fact, we will show that m has to be an odd integer in order to avoid 
collisions between masses traveling on different loops. By symmetry, each 





flfc sin kt + bk cos kt , 



k 
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pair of loops intersects at two diametrically opposite points. If m is even, 
then the Lagrange symmetry implies that each mass has a "partner" at the 
diametrically opposite point on the same loop. Hence, at the moment when 
one of the masses crosses an intersection with another loop, its partner crosses 
the opposite intersection. However, the masses on different loops are related 
by the rotations (jH]), and so at the same moment this other loop also has a 
mass at each of these two intersections, and there will be two simultaneous 
collisions (and, by symmetry, collisions for other pairs of loops as well). 

By the same reasoning, if there is an odd number m of masses on each 
loop, at the moment a mass on one loop crosses an intersection with another 
loop, the corresponding mass on this other loop crosses the other intersection 
between the two loops. Thus no two masses cross an intersection at the same 
time, and no collisions occur. 

Table [1] gives Fourier coefficients found by our action minimization for 
m = 1,3,5,7, normalized with ai = 1. As for m = 1, we checked that the 
force and acceleration agreed to within 5 decimal figures at all times. Movies 
of these orbits are available at [TB] ; a movie of the case m = 1 also appears 
at 0. 

The special case when m is a multiple of 3 has some additional symme- 
tries. By Lagrange symmetry, the set of masses is fixed under translating 
time by 2tt/3. Therefore, at all times the set of masses is now fixed not 
just under the four transformations (jS]), but also under the cyclic symmetry 
X y ^ z ^ X. These transformations generate a 12-element symmetry 
group, consisting of half the orientation-preserving rotations of the cube (for- 
mally, a subgroup which is isomorphic to the alternating group A4). In terms 
of linear transformations, these are the matrices 



(10) 



with the restriction that the number of —Is in each matrix is even. 

This symmetry has amusing consequences. Recall that the moment of 
inertia tensor I and the quadrupole moment Q are defined as 

J = ^mi|xip(l - TTxJ, Q = ^mi|xi|^(37rx, - 1) 

i i 

where tTx denotes the projection operator onto the subspace parallel to x, 
and 1 denotes the identity matrix. Then / and Q must commute with each 
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k 


Ofc (m = 1) 


Ofc (m = 3) 


Ofc (m = 5) 


ttk {m = 7) 


1 


1.00000 


1.00000 


1.00000 


1.00000 


3 


.03282 


-.04629 


-.03991 


-.03335 


5 


-.00098 


-.00472 


.00359 


.00885 


7 


-.00036 


-.00269 


.00161 


.00445 


9 


-.00003 


.00056 


-.00168 


-.00334 


11 




.00010 


-.00043 


-.00039 


13 




.00007 


-.00028 


-.00028 


15 




-.00002 


.00012 


-.00013 


17 






.00001 


-.00011 


19 






.00001 


-.00008 


21 








-.00005 



Table 1: Fourier coefficients ak for the cubic orbits with ai normahzed to 
unity. There are m masses on each loop for a total of 4m masses. 

matrix in the symmetry group ( fTOl) . However, since ( fTOl) is an irreducible 
representation of ^4 [7], the only such matrices are scalars, i.e., multiples of 
the identity. Moreover, since Q has zero trace, it must be identically zero. 
Thus, whenever m is a multiple of 3, / and Q are just as they would be for 
a sphere, indicating that the mass distribution of these orbits is spherical to 
second order. (The moment of inertia is varies with time, however, so these 
orbits are not counterexamples to Saari's conjecture [19] that the only orbits 
with constant moment of inertia are central configurations.) 

Finally, we note an interesting property of these orbits. As m grows large, 
the Fourier coefficients a^, normalized with ai = 1, appear to converge (albeit 
slowly) to asymptotic values for each k: Table [2] gives the ffist 5 nonvanishing 
Fourier coefficients for m = 61,81, 101, and 121. Based on this observation, 
we conjecture that as m — > 00, these orbits approach a fixed asymptotic 
shape. A similar situation seems to hold for the figure-8 [20l [17] . 

3 The Criss-cross 

The three-dimensional periodic orbits we discuss above are unstable. In 
contrast, we present now a two-dimensional periodic orbit for three masses 
which we call the "criss-cross." This orbit was ffist obtained in 1976 by 
Henon [9] as a member of a family of finite angular momentum orbits for 
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k 


(m = 61) 


Ofc (m = 81) 


ttk {m = 101) 


O'k {fn = 121) 


1 


1.00000 


1.00000 


1.00000 


1.00000 


3 


-.02302 


-.02285 


-.02265 


-.02235 


5 


.01054 


.01025 


.01008 


.00991 


7 


.00544 


.00522 


.00507 


.00494 


9 


-.00181 


-.00175 


-.00167 


-.00156 



Table 2: Fourier coefficients Ofc for the orbits with cubic symmetry with ai 
normahzed to one. Here n = 4m, with m masses on each loop. 

equal masses, which are an extension of the zero angular momentum Schubart 
orbit. In 1993 it was rediscovered by Moore by searching for orbits with a 
particular braid type [15], and later shown numerically to exist for different 
masses by Nauenberg [17]. More recently a rigorous proof of existence has 
been given by Chen [1] as part of a family of retrograde orbits. 

Henon [9] studied the linear stability of the extended Schubart family, and 
predicted that it contains a finite interval of periodic orbits which are stable 
in three dimensions. The extent of this interval, and whether it includes the 
criss-cross orbit, was not calculated. However, even if an orbit is stable to 
ffist order, it is entirely possible that is not dynamically stable. Here we 
provide some numerical evidence that the criss-cross orbit is in fact stable to 
perturbations, not just in the plane, but in three dimensions. 

The Fourier expansion for the coordinates is given by 

Xi(t) = ^ Qi^k COS kt, yi = ^ bi^k sin kt (11) 

for i = 1,2,3. Under the transformation t — > t + vr, these coordinates 
change sign, and therefore the sum is restricted to odd values of k. For 
equal masses, the Fourier coefficients also satisfy the symmetry condition 
Ci2,k = &2,fc = ±0,1, fc and 63, ^ = ±03,^, where the plus sign applies 

for k = 3,7, 11, 15 . . . and the minus sign applies for k = 1, 5, 9, 13 . . . Here 
masses 1 and 2 are those orbiting each other in the center, while mass 3 
orbits them both in the retrograde direction. 

To start the Fourier gradient search of the minimum of the action integral, 
we take for the initial Fourier coefficients ai^i = 1, 61,1 = 0, and 03 1 = — 1. 
The resulting coefficients, which provide Fourier sums that converge to 5 
decimal figures, are shown in Table 3. In Fig. [3] on the left, we show the 
orbits for three equal masses at 16 equal time intervals. One of the three is 
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7 

k 




bi,k 




1 


1.09764 


0.10896 


-0.98868 


3 


-0.02809 


0.03251 


-0.00442 


5 


0.00724 


-0.00376 


-0.01100 


7 


-0.00121 


0.00131 


-0.00010 


9 


0.00040 


-0.00029 


-0.00069 


11 


-0.00010 


0.00010 


-0.00001 


13 


0.00003 


-0.00003 


-0.00006 


15 


-0.00001 


0.00001 


0.00000 



Table 3: Fourier coefficients for the equal mass criss-cross orbit. 

nearly circular and rotates in the opposite direction of the other two orbits. 
The initial conditions for the three masses are 

xi = (1.07590,0), xi = (0,0.19509) 
X2 = (-0.07095,0), ±2 = (0,-1.23187) 
X3 = (-1.00496,0), ±3 = (0,1.03678) 

Similar orbits exist for other mass ratios as well. On the right of Figure [3] we 
show the corresponding orbit with three masses in the ratio 1:2:3, showing 
that similar orbits exists for mass ratios far from 1. 

In Fig. Hlon the left, we show the equal mass criss-cross orbit during 40 
periods for an initial small deviation 6x = .001 of the initial position of one 
of the masses, indicating that the resulting orbit is quasi-periodic. For longer 
times the orbit fills out the circle. On the right side we show a corresponding 
quasi-periodic orbit which results for a somewhat larger initial deviation 
Sx = .005 of the initial position. Similarly, a deviation Sz = .005 along 
the 2;-axis results in the beautiful quasi-periodic orbit in three dimensions 
shown edge-wise in Fig. El To exhibit this behavior we have magnified the 
axis normal to the plane of the original periodic orbit by a factor of 50. 

These sample results suggest that this orbit is dynamically stable with 
respect to a fairly wide range of perturbations, both in the plane and in three 
dimensions. This stability, and the fact that similar orbits exist for a wide 
range of mass ratios, suggests that triple systems of this type may actually 
exist in astronomical systems. 
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Figure 3: On the left, the criss-cross orbit for three equal masses shown at 
16 equal time intervals. Starting with the masses aligned horizontally, at the 
second time interval they lie at the vertices of an isosceles triangle. On the 
right, the corresponding criss-cross orbit with three masses in the ratio 1:2:3. 




Figure 4: On the left, the equal mass criss-cross orbit during 40 periods, 
for an initial small deviation 6x = .001 of the initial position of one of the 
masses. On the right, the same orbit with a larger perturbation 6x = .005 
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Figure 5: Perturbing the criss-cross along the z direction leads to this beau- 
tiful quasiperiodic orbit, indicating that the criss-cross is dynamically stable 
in three dimensions. The vertical component of motion has been exaggerated 
by a factor of 50 for purposes of illustration. 
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